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ABSTRACT 

The gravitational- wave and accretion driven evolution of the angular veloc- 
ity, core temperature, and (small) amplitude of an r-mode of neutron stars in 
low mass X-ray binaries and similar systems is investigated. The conditions re- 
quired for evolution to a stable equilibrium state (with gravitational wave flux 
proportional to average X-ray flux) are determined. In keeping with conclusions 
derived from observations of neutron star cooling, the core neutrons are taken to 
be normal while the core protons and hyperons and the crust neutrons are taken 
to be singlet superfluids. The dominant sources of damping are then hyperon 
bulk viscosity (if much of the core is at least 2-3 times nuclear density) and (e-e 
and n-n) shear (and possibly magnetic) viscosity within the core-crust boundary 
layer. It is found that a stable equilibrium state can be reached if the superfluid 
transition temperature of the hyperons is sufficiently small (< 2 x 10 9 K), al- 
lowing the gravitational radiation from Sco X-l and several other neutron stars 
in low-mass X-ray binaries to be potentially detectable by the second generation 
LIGO (and VIRGO) arrays. 

Subject headings: accretion, accretion disks - - dense matter - - gravitational 
waves — stars: neutron — X-rays: binaries 



1. Introduction 

In this Letter, we present the results of an investigation of the evolution of rapidly 
rotating accreting neutron stars under the influence of their emission of gravitational ra- 
diation. We employ a modification (Wagoner, Hennawi Sz Liu 2001) and extension of the 
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two- component (equilibrium plus perturbation) model of the star introduced by Owen et al. 
(1998) and also employed by Levin (1999), but restrict the analysis to small perturbations. 
These are assumed to be in the form of r-modes (Andersson 1998; Friedman & Morsink 
1998; Lindblom, Owen & Morsink 1998; Andersson, Kokkotas & Schutz 1999), which radi- 
ate mainly via Coriolis-driven velocity perturbations rather than the density perturbations 
of the less powerful f-modes. We must allow for large uncertainties in many of the relevant 
properties of neutron stars, such as the superfluid transition temperatures and the properties 
of the core-crust boundary layer. Details will be presented in a subsequent paper. 

After developing a general formalism, we shall focus on conditions in which the neu- 
tron star angular velocity (and thus gravitational wave frequency) evolves slowly toward an 
equilibrium state, in which the rate of accretion of angular momentum from the surrounding 
disk is balanced by its rate of loss via gravitational radiation. If this equilibrium is achieved, 
the observed flux of gravitational radiation can be shown to be proportional to the observed 
flux of X-rays from the accretion (Wagoner 1984; Bildsten 1998). 

One of our longer term goals is the development of parameterized expressions describing 
possible time evolutions of the gravitational- wave frequency and amplitude, to facilitate 
detection by LIGO, VIRGO, and similar laser interferometer detectors. The brightest low 
mass X-ray binaries (LMXBs) thought to contain a neutron star are the prime targets. 

2. Dynamical and Thermal Evolution 

In this exploratory investigation, it will be sufficient to consider a Newtonian neutron 
star in equilibrium (with equatorial radius R) which is perturbed by a nonaxisymmetric 
infinitesimal fluid displacement £ = f{r,6)e l{ - m<l>+ ' 7t > ~ aR, with a « 1. Based on the 
work of Friedman & Schutz (1978a) and Levin & Ushomirsky (2001a), the total angular 
momentum J of the star can be decomposed into its equilibrium angular momentum J* and 
a perturbation proportional to the canonical angular momentum J c . That is, 

J = J*(M, n) + (1 - Kj)J c , J c = -K c a 2 J* , (1) 

where M is the mass and fl is the (uniform) angular velocity of the equilibrium star. All 
constants K( ) will be dimensionless, with Kj ~ K c ~ 1. 

In classical mechanics, the action I = E/uo of any normal mode of a set of oscillators 
(with frequency a;) is an adiabatic invariant. For a fluid, the analogous quantity should 
be E c /cu, where E c is the canonical energy of the perturbation in the corotating frame and 
uj = a + mQ is its frequency in that frame. However, we also have the general relation 
E c = -(uj/m)J c (Friedman & Schutz 1978a). Therefore, following Ho & Lai (2000), we 
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assume that the canonical angular momentum is also an adiabatic invariant, and should 
therefore be unaffected by the slow rate of mass accretion. Thus it obeys the usual relation 
(Friedman & Schutz f978b) 

dJ c /dt = 2J c [(F g (M,n)-F v (M,n,T v )] , (2) 

where F g is the gravitational radiation growth rate and F v is the viscous damping rate. The 
latter also depends upon a spatially averaged temperature T v (t). 

On the other hand, conservation of total angular momentum requires that 

dJ/dt = 2J c F g + j a (t) , (3) 

where J a = j a M is the rate of accretion of angular momentum. The mass is accreted with 
specific angular momentum j a at a rate M{t). 

Combining these equations then gives the dynamical evolution relations 



~^ = F g -F v + [K j F g + (l-K j )F v ]K c a 2 -(jf)M(t), (4) 



(£)§ = -2[K J F g + (l-K j )F v }K c a 2 + 



(jo - 3*) 

J * 



M(t) ; (5) 



where /*(M, fi) = dJ^/dVt and j*(M, Q) = dJ*/dM. In keeping with the fact that it is 
sufficient to also work to lowest order in fi/fi(max) (as well as the relativity parameter 
GM/Rc 2 ) in this exploratory investigation, we take j a — j* = j a and J* = i*f2. We also note 
that K c = 0.094 (Owen et al. 1998). (The value of Kj is unimportant, since we will see that 
F v remains very close to F g .) 

Finally, thermal energy conservation for the entire star gives our third evolution equation 

f BT dT 
J —c v dV = C(T) — - 2E C F V (T V ) + K n (M)c 2 - L U (T U ) , (6) 

where the corotating frame canonical energy E c = K e QJ*a 2 , with K e = K c /3. Since all 
constituents (electrons, nucleons,. . .) are degenerate, the specific heat at constant volume (c v ) 
is essentially the same as that at constant pressure. The terms on the right hand side of this 
equation represent viscous heating, pycnonuclear reactions and neutron emissions in the inner 
crust (proportional to a time-averaged mass accretion rate), and neutrino luminosity. The 
hydrogen/helium burning rate is assumed to be balanced by the surface emission of photons 
(Schatz et al. 1999), especially at the large accretion rate M = 1O~ 8 M yr" 1 (roughly 1/3 
the Eddington rate, appropriate to our primary targets) that we shall adopt. The mass 
accretion rate can be estimated from accretion energy conservation. The photon luminosity 
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arising directly from the accretion is L acc ps (GM/R)M(t), for a slowly rotating neutron star 
with a negligible magnetosphere. 

From now on we shall take the perturbation to be due to the dominant / = m = 2 
r-mode, in which case the gravitational wave frequency f gw = (4/3) / ns = 2fl/3n. In order to 
facilitate comparison with previous results, we shall adopt the neutron star model of Owen 
et al. (1998) (p oc p 2 , M = 1.4M , and Q c = (ttG(p)) 1/2 = 8.1 x 10 3 rad/s) for numerical 
work. Then the gravitational radiation growth rate of this mode is 

F g = fi 6 /v » V = 3 - 26 s > ^ = Q / Q c , (7) 
which should also dominate that due to the f-modes. 

Comparison of observations of thermal emission from isolated neutron stars with com- 
puted cooling histories (Tsuruta et al. 2002) has led Kaminker, Yakovlev & Gnedin (2002) to 
conclusions regarding the most likely state of the constituents, which we tentatively adopt. 
Specifically, the maximum values of the (density dependent) superfluid transition temper- 
atures are taken to be (a) T n < 10 8 K for the (triplet) core neutrons (so they will remain 
normal), (b) T p > 5 x 10 9 K for the (singlet) core protons, and (c) T n > 5 x 10 9 K for 
the (singlet) inner crust neutrons. In what follows we shall also assume that the thermal 
conductivity timescales r th are short enough to give relations T V (T) and T„{T) between these 
three spatially averaged temperatures that appear in equation (6). Typically, r th ~ 0.1 — 1 
year in the core and r t h ~ 10 — 10 2 years in the inner crust (Brown 2000; Baiko, Haensel & 
Yakovlev 2001). 

Then in the temperature range of interest (10 8 < T < 10 9 K), the viscous damping rate 
of this mode is 

F v = F sh (T) + F bl (Q,T) + F hb (Q,T) . (8) 

Even if the neutrons were not normal, the contribution to this damping rate from the mutual 
friction between a neutron superfluid and the superconducting proton — relativistic electron 
fluid (Lindblom & Mendell 2000) would be negligible. 

The first term is produced by the ordinary shear viscosity throughout the core. Adding 
the contribution of the n-n scattering (Owen et al. 1998) to that of the e-e scattering (Lind- 
blom & Mendell 2000) gives 

F sh = l/(r sh Ti) , T 8h = 0.69 x 10 6 s , T 8 = T/10 8 K . (9) 

The second term is produced by the ordinary shear and magnetic viscosity in the crust- 
core boundary layer. For a normal fluid, this damping rate has been calculated by Andersson 
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et al. (2000), Bildsten & Ushomirsky (2000), Lindblom, Owen & Ushomirsky (2000), Rieu- 
tord (2000), and Wu, Matzner & Arras (2001) neglecting magnetic effects, while Bildsten 
& Ushomirsky (2000), Mendell (2001), and Kinney & Mendell (2002) included them. For 
superffuid neutrons and protons, this rate was calculated by Kinney & Mendell (2002). We 
shall approximate the damping rate for our case (normal neutrons and superffuid protons in 
the boundary layer) by adding the rate of Kinney & Mendell (2002) for a normal uncharged 
fluid to their superffuid rate when only the proton vortices (magnetic flux tubes) are pinned 
to the crust. The magnetic damping is due to the interaction of the field with the effective 
current density of the proton vortices, producing cyclotron-vortex waves. Although there is 
no mutual friction because no neutron vortices are present, we also neglect the small cou- 
pling between these fluids produced by n-e scattering. We then obtain (after an approximate 
algebraic simplification) 



where B (Gauss) is the radial magnetic field. The slippage factor S s is the fractional degree 
of pinning of the vortices in the crust (Kinney & Mendell 2002); and S 2 S = (2S 2 + S 2 )/3, 
with the slippage factor S n the fractional difference in velocity of the normal fluid between 
the crust and the core (Levin & Ushomirsky 2001b). Although both slippage factors were 
defined to be at most unity, we also let them contain the uncertainties in our model. [If 
the pinning were dominated by neutron vortices, the coefficients in equation (10) would be 
(2 — 6) x 10 3 times larger, resulting in no growth of the mode unless «S S <C 1.] In obtaining 
most of the numerical results below, we take the magnetic contribution in equation (10) to 
be negligible, which requires that B <6x 10 10 [1 + 2(S n /S s ) 2 )] 2 nT 8 2 G. 

The third term arises from the bulk viscosity produced by out-of-equilibrium hyperon 
reactions (which dominate that produced by direct and modified Urea reactions). This has 
been studied by Jones (2001a,b) and Lindblom & Owen (2002) for normal nuclear matter and 
by Haensel, Levenfish & Yakovlev (2002) for superfluids. We employ the results of Lindblom 
& Owen (2002) for n + n ^ p + £ _ . However, with our assumptions about superfluidity, the 
reaction whose rate is least reduced by superffuid phase space blocking (but is more difficult 
to calculate) should be n + n ^ n + A. These hyperons should be present at densities 
p > (5 — 8) x 10 14 g cm -3 , which is achieved over a large fraction of the core of relevant mass 
[M = (1.3 — 1.6)M©] neutron star models for many nuclear equations of state (which are 
softened by their presence) (Balberg, Lichtenstadt & Cook 1999). Employing the superffuid 
reduction factor R hb (T/T h ) of Haensel, Levenfish & Yakovlev (2002) for normal nucleons and 
a hyperon with singlet superffuid transition temperature Th, we obtain 



F bl = 3.15 x 10~ A Sl s Vn/T 8 + 3.71 x 10~ 8 S 2 y/B s 



-i 



(10) 



F hb = f hb tfr{T)/[l + (2fir(T)/3) 2 ] , r(T) = t^ 2 / R hb {T/T h ) , 



(11) 
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where 1/r is the reaction relaxation rate. The constants t « 0.15 s and t\ ps 0.01 s are 
obtained by fitting the results of Lindblom & Owen (2002). The factor f hb allows for the 
difference in the rate of the two reactions and the fraction of the neutron star that has a 
density high enough to allow the reactions. Since the dominant uncertainty in equation(ll) 
is due to that of Th, we shall set fhb — 1, implying that the neutron star is sufficiently massive 
to have a central density well above the indicated hyperon threshold. 

The temperature dependences of these contributions to the viscous damping rate are 
shown in Figure 1. Now that we have specified all properties in the equations (4) and (5) 
of evolution of a(t) and Q(t), we can consider the thermal evolution [equation (6)]. With 
the neutrons normal, they dominate the specific heat, giving C(T) = 1.5 x 10 38 T 8 erg/K. 
(The electron contribution is about 15 times less.) We also take the nuclear heating constant 
K n = 1 x 10~ 3 (Brown 2000). 

The neutrino luminosity is taken to be 

L u = L du T% 'R du (T ] /T p ) + L mu Tg R mu (T /T p ) + L ei T% + L nri T 8 8 + L cp Tj . (12) 

The proton superfluid reduction factors for the direct Urea reactions (Rdu) and the modified 
Urea reactions (R mu ), as well as most of the following constants, are obtained from the 
review of Yakovlev, Levenfish & Shibanov (1999). The other terms represent (inner crust) 
electron-ion and (core) neutron-neutron neutrino bremsstrahlung, and Cooper pairing of 
(inner crust) neutrons. We take L du = f du x 10 8 L mu , where f du includes the fraction of the 
neutron star that is above the density threshold for the direct Urea reactions (comparable to 
the above threshold for the hyperon reactions). The constants L mu 1.0 x 10 32 erg/s and 
L e { ps 9.1 x 10 29 erg/ s are obtained by fitting the results of Brown (2000) (who did not consider 
the other processes) for normal and superfluid neutron stars. Finally, L nn ps 0.01L mu , while 
L cp ps 9 x 10 31 erg/s is proportional to the length scale of variation of the neutron superfluid 
transition temperature within the inner crust (taken to be 100 m) (Yakovlev, Kaminker & 
Gnedin 2001). 

We are interested in the evolution of neutron stars after they have been spun up to the 
point where the gravitational radiation growth rate has become equal to the viscous damping 
rate: 

F g (n ,M ) = F v (n ,M ,T ) = F . (13) 

This equality defines our initial state. Before that time, we see from equation(2) that any 
intrinsic perturbation could not grow from its (infinitesimal) value a min . The initial tem- 
perature T is then determined by the vanishing of equation (6), with the nuclear heating in 
the inner crust balanced by the neutrino emission (Brown 2000). This temperature is most 
sensitive to the direct Urea factor f du , but only varies by 14% over the range < f du < 1. 
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Consistent with our assumption that the hyperon reactions are operative, we shall adopt 
fdu = 0.1. 

For a typical value f2 ~ 0.3, the time scale r = l/F ~ 10 4 sec. Two other key time 
scales are that due to cooling and accretion, 

_ = F = L ^ „ _J_ L = F = (k\lM)~^— (14) 
r c ~ c "C(To)T 103 yr' r a " * a ~ { M} W yr " U4J 

Since we will be concerned with time scales At <C M / (M) > 10 8 yr, we can take M(t) = M . 

In contrast to the initial state, the equilibrium state X\ of our dynamical variables 
X l (t) = {a,Q,T} is defined by the vanishing of the evolution equation (3), in addition to 
equations (2) and (6). From equation (4) or (5), one can see that the equilibrium amplitude 
is given by 

1/2 

~ (10" 6 - 10" 5 ) . (15) 



a P = 



F a 



_2K c F g {tt e )_ 

These values of a e are much less than the saturation amplitude of the r-mode instability 
(Arras et al. 2002). 

The linearized analysis of Wagoner, Hennawi & Liu (2001) shows that stability of the 
equilibrium requires that 

1 dL v \ ( 1 dR, 



> 7^ > , (16) 



L v dT J e \F V dT 

1 /2 

assuming that \d/dT\ ~ 1/T. If so, oscillations of frequency f a ~ K r cn e F e are damped at a 
rate fa ~ K r a 2 e F ei where F e is the equilibrium value of F v = F g and K r = 2K e Q e J /C e T e ~ 
10 5 is the ratio of rotational to thermal energy. We have also used the fact that the viscous 
heating term in equation (6) is typically at least ten times greater than the nuclear heating 
term. In practice, the first inequality in equation (16) is always satisfied, if the key require- 
ment dF v /dT > holds (which is the case if F v is dominated by the hyperon bulk viscosity, 
as shown in Figure 1). This requirement applied to the initial state also guarantees slow 
evolution, as we shall see below. 



3. Conclusions 

In Figure 2 we show the critical curve fl(T 8 ) given by F g = F v , for three choices of the 
key parameters T h ,S ns . Also shown are the initial state fl ,T /10 8 K and the equilibrium 
state fi e ,T e /10 8 K, with T e > T = 3.32 x 10 8 K. For T k >3x 10 9 K, these results are 
independent of T h . For T h < 1 x 10 9 K, these results are independent of the slippage factors 
S n and S s (and Q e does not increase greatly as T h decreases in this range). 
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It can be shown that the evolution from the initial state is also controlled by the sign 
of (dF v /dT) , which is equal to the sign of the slope of the critical curve. If the slope is 
negative [case (a)], there will be a thermal runaway with a growth rate f Q (given above) 
that is of the same magnitude as found by Levin (1999). If the slope is close to zero [case 
(b)], there will initially be overstable oscillations of the type found by Wagoner, Hennawi & 
Liu (2001). If the slope is positive [case (c)], the oscillations of the growing amplitude are 
damped out on a timescale r c (as shown in Figure 3), after which it slowly increases to its 
equilibrium value [a e = 1.9 x 10 -6 for the parameters of case (c)], along with f2 and T. The 
time required to reach equilibrium is At (AQ/Q)r a . Throughout, F g remains very close 
to F v , so the evolution is along the critical curve. 

The maximum rotation rate (due to shedding) of our chosen neutron star model is 
fmax — (2/3)(fi c /27r) = 856 Hz. Coherent oscillations in Type 1 X-ray bursts have been 
observed at frequencies F < 590 Hz (van der Klis 2000), which would then correspond to 
Q < 0.46 if they represented the rotation rate. There is evidence that some of these are the 
first harmonic, in which case the highest spin frequency is 350 Hz (van der Klis 2000). (For 
isolated neutron stars, / < 642 Hz.) 

Observations of the luminosity of LMXB's in their quiescent (low accretion rate) phase 
constrain the r-mode amplitude if they are in the equilibrium state we have considered 
(Brown & Ushomirsky 2000), where the r-mode viscous heating balances the neutrino energy 
loss. This will also be addressed in a subsequent investigation. 

Our main conclusion is that evolution to a stable equilibrium state can occur if (a) a 
significant fraction of the core of the neutron star is above the threshold for hyperons, (b) 
their superfluid transition temperature T/,<2x 10 9 K, (c) the core neutrons near the crust 
are not a superfluid whose vortices are strongly pinned to the crust, and (d) the magnetic field 
is not too strong in that core-crust boundary layer. If Sco X-l has been spun up by accretion 
to such a stable equilibrium state, it should be detectable by the second generation LIGO 
detectors. (However, its spin period remains unknown.) When signal recycling ('narrow- 
banding') is employed, a few additional LMXB's may also be detectable (Cutler & Thorne 
2002). 

This work was supported in part by NSF grant PHY-0070935. We benefitted from 
many discussions during the 2000 program on Spin and Magnetism in Young Neutron Stars 
at the Institute for Theoretical Physics, U.C. Santa Barbara. Joseph Hennawi and Jingsong 
Liu provided important input to an early stage of this investigation. Thanks also goes to 
Lee Lindblom and Greg Mendell for very helpful discussions. 
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Fig. 1. — The dependence on temperature (Tg) of the three contributions to the damping 
rate F v : core shear viscosity (long dashed), boundary layer viscosity for B < 10 9 Gauss 
(short dashed) and B = 10 11 Gauss (short and long dashed), and hyperon bulk viscosity 
(solid). The model chosen has T h = 2 x 10 9 K, S n = S s = 0.2, and Q = 0.30. 
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Fig. 2. — The relation bewtween f2 = f2/f2 c and T 8 on the critical curve defined by F t 
is shown for the choices (a) T h = 3 x 10 9 K, S ns = 1.0 (long dashed), (b) T h = 3 x 10 1 
K, 5 ns = 0.1 (short dashed), and (c) T h = 1 x 10 9 K, S ns < 0.3 (solid). The magnetic 
boundary-layer viscosity is assumed negligible. Also shown are the initial and equilibrium 
states (at the higher temperature on each curve). 
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Fig. 3. — The early evolution of the r-mode amplitude a, for case (c) in Figure 2. The initial 
amplitude was chosen to be oio = 1CT 12 , and t* — 1 year. 



